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Abstract 

We present non-perturbative results for the constants needed for on-shell 0{a) improvement 
of biUnear operators composed of Wilson fermions. We work at /3 = 6.0 and 6.2 in the quenched 
approximation. The calculation is done by imposing axial and vector Ward identities on correlators 
similar to those used in standard hadron mass calculations. A crucial feature of the calculation is 
the use of non-degenerate quarks. We also obtain results for the constants needed for off-shell 0{a) 
improvement of bilinears, and for the scale and scheme independent renormalization constants, Za, 
Zy and Zs/Zp. Several of the constants are determined using a variety of different Ward identities, 
and we compare their relative efficacies. In this way, we find a method for calculating cy that gives 
smaller errors than that used previously. Wherever possible, we compare our results with those of 
the ALPHA collaboration (who use the Schrodinger functional) and with 1-loop tadpole-improved 
perturbation theory. 
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I. INTRODUCTION 

Symanzik's improvement program is a systematic method for reducing discretization 
errors in lattice simulations ||l], |^ . One must improve both the action and external operators 
by the addition of appropriate higher dimension localized operators. Complete removal of 
discretization errors at a given order in the lattice spacing, a, requires a non-perturbative 
determination of the coefficients (the "improvement constants") of the higher dimension 
operators. A key ingredient in the practical implementation of the improvement program is 
the development of methods for such non-perturbative determinations. 

The ALPHA collaboration has exploited the connection between 0{a) discretization er- 
rors and chiral symmetry to develop non-perturbative methods for the calculation of some of 
the 0{a) improvement constants (those for the action and some of the local fermion bilinear 
operators) [^ ^, ^, ^ . Their approach is based on the imposition of axial and vector Ward 
identities. It also determines the renormalization-scale independent normalization constants 
Z^, Zy and Z^/Zp, as originally observed in Ref. ^ This non-perturbative determination 
of improvement and normalization constants is of considerable practical importance, as un- 
certainties in these constants can be a significant source of error in lattice calculations of 
matrix elements. 

In Ref. ^ we showed how to extend the method of the ALPHA collaboration to determine 
all the 0{a) improvement constants for bilinears.^ The extension involves the enforcement 
of Ward identities for massive, non-degenerate quarks, rather than in the chiral limit, and is 
a generalization of the method of Ref. |lT|. Results of a pilot simulation at /5 = 6 (quenched) 
suggested that the method was practical. This simulation had the drawback, however, 
that it was done using tadpole-improved, rather than non-perturbatively improved, Wilson 
fermions. Thus a clean separation of sources of error was not possible. 

In this paper we present results of a more extensive investigation of the method. We 
use the non-perturbatively improved action, taking the non-perturbative value for the 



Sheikholeslami-Wohlert (or "clover" ) coefficient csw |T2[ from the work of the ALPHA col- 
laboration 0. Thus the errors after improvement should be of Ola'^). We study the scaling 
behavior of improvement and normalization constants by carrying out the calculation at two 
values of the lattice spacing, (3 = 6 and 6.2 (quenched). We also extend previous work by 
determining the improvement coefficients for the operators which vanish by the equations 
of motion ( "equation-of-motion operators"). These contribute only to off-shell matrix ele- 
ments, and thus are not of direct physical relevance, but they do contribute to the Ward 
Identities at non-zero quark masses. 

As already noted, several of the improvement and renormalization constants that we de- 
termine have been previously obtained by the ALPHA collaboration. An important differ- 
ence in the implementation of the improvement conditions is that the ALPHA collaboration 
uses Schrodinger functional boundary conditions with sources on the boundary, while we 
use periodic boundary conditions with standard sources for quark propagators designed to 
improve overlap of local operators with hadronic ground states. This means that the results 
for improvement constants will differ at 0{a) and the normalization constants will differ at 
0(0^). One of the aims of our study is to compare results from the two approaches, since 
this gives an indication of the importance of the neglected higher order terms. We can also 



^ Other approaches that allow one to determine, in principle, all the improvement and normalization con- 
stants have been suggested in Refs. ^, |l^. 



get some idea of the relative effectiveness of the two approaches. 

The organization of this paper is as follows. In the following section we briefly reca- 
pitulate the theoretical background to our method, and give a general description of our 
implementation. Sec. |T| contains a summary of our simulation parameters. In Sec. |V|, we 
present our final results, and discuss their implications. We reserve a detailed discussion of 
the calculation of the individual improvement coefficients for Sees. |V} |X-14 We close with 
some conclusions in Sec. pClll|. Three appendices collect the tadpole- improved perturba- 



tive results which we use for comparison with our non-perturbative estimates, the tree-level 
definitions of the improvement constants, and a discussion of exceptional configurations. 

II. WARD IDENTITIES: THEORETICAL BACKGROUND 

On-shell improvement of bilinear operators at 0{a) requires both the addition of extra 
operators, 

{Ai)f, = A^ + acAd^P 
{Vi)fM = V^ + acvdyT^y 

(Tj)^, = T^, + acTid^V, - d,V^) (1) 

Pi = P 
Si = S, 

and the introduction of the following mass dependence 

Og^) ^ Z'ai^ + boam,,)Of'\ (2) 

= Z%{l + hoam,,)Of\ (3) 

Here {ij) (with i ^ j) specifies the fiavor, and O = A,V, P, S, T. The Zq are renormalization 
constants in the chiral limit, rriij = {rrii + mj)/2 is the average bare quark mass,^ and rhij 
is the quark mass defined in Eq. (|T5D using the axial Ward identity (AWI). There are 
yet other improvement constants needed in order to extend the analyses to fiavor-neutral 
bilinears (i = j) and to full QCD. These extensions are discussed in Ref. |1^, but are not 
relevant here. Note that, except in Appendix ^, we have set the Wilson parameter r equal 
to unity. 

When improving the theory to 0{a), one still has freedom in defining the cq and the ho. 
For example, in general, they can depend on the correlators used to define them and on the 
quark mass. We shall consistently use the value in the chiral limit as it is the simplest choice 
and is also the one made in previous work by other collaborations. The correlators used to 
define them are discussed in subsequent sections. 

To avoid confusion, we stress that the coefficients ho differ from the ho used by earlier 
authors. In particular, at the level of 0{a) improvement, one has 
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[Z\Z%lZl)ho . (4) 



^ The bare quark masses are ami — l/2Ki — 1/2kc, k being the hopping parameter in the Sheikholeslami- 
Wohlcrt action and Kr its value in the chiral limit. 



The analogous relation between m and m is given in Eq. 

Improvement can be achieved by imposing the generic axial Ward identity 

iSS^''^ O'Sffiy) J^'^HO)) = {60f}g{y) J(3^)(0)) (5) 

for enough choices of J, O, and y to determine all the relevant improvement and scale 
independent normalization constants. This should then guarantee that the identity holds 
up to corrections of 0{a?') for other choices of J and y. Here 50 is the bilinear which results 
from the axial variation of O in the continuum {A^ ^-> V^, S ^^ P, and T^,y -^ ^fj.upaTpa), 
and the variation of the action under an axial rotation is 



^^(12) ^ ^(12) / ^4^ 



(12) _fi (A^ \{12) 



i2mu)iPi,offr'^-dMi,off 



(6) 



The point y lies within the domain, V, of the chiral rotation, while the source J is located 
outside V. 

To implement Eq. (^) away from the chiral limit, it is not sufficient to use the on-shell 
improved bilinears. Or, defined in Eqs. (|,^. One must also include dimension 4 operators 
which vanish by the equations of motion, and this has been anticipated in the use of the 
subscript off. As noted in Ref. ^ there is one such operator with the appropriate symmetries 
for each bilinear: 

0^;% = Of^ - a\c'oE^S^ , (8) 

^g^') = ^^'^rvi)^(^') - ^^"^VvT^^^^ . (9) 

In the equation-of-motion operators Eq, F is the Dirac matrix defining O, and Wipj = 
{p + 'mj)ipj + Ola"^) is defined to be the full 0{a) improved Dirac operator for quark flavor 
j (see Appendix ^. This ensures that Eo gives rise only to contact terms, and thus cannot 
change the overall normalization Zq. The factors multiplying Eq are chosen such that, at 
tree level, Cq = 1 for all Dirac structures as shown explicitly in Appendix ^. 

For practical applications, it is useful to express the Ward identity in terms of on-shell 
improved operators. The equation-of-motion operators contribute only when the operators 
P (contained in SS) and O, in the l.h.s. of Eq. (||), coincide. The 75 in Pi, off changes Oj^off 
to SOj^off, and so, up to O^a^) corrections, these contact terms are proportional to the r.h.s. 
of Eq. (^). After rearrangement, one finds 

iU'.mor>(y,,^J'^-^m_ ZS? ^„£k±^-^^^o(a=) (10) 



where 



{60\''\y„y) J(^^)m Zf^zf) 



6Sj{x) ^ 2muPi''\x) - d,{AjY^'\x) . (11) 



This is the form of the AWI which we enforce {i.e. for some choice of J we fit to a range in 
y, neglecting O(a^) contributions) in order to determine the improvement constants. Note 
that the mass multiplying the c' coefficients is m and not m. Also, for brevity, mention of 
the O(a^) terms in all equations is henceforth omitted. 



To highlight the dependence on quark masses, the r.h.s. of Eq. ([T0|) can be written as 
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^■0 r _ ^ ^ ^ -| c' + C' 

1 + absomi^ - abAm^ - abom23 + a ^ 



and, in the special case mi = m2 relevant to our numerical study, as 



+ a-^^^^12 , (12) 
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ami • (13) 



Here we have defined rhi = 'mij\m=m,, ^-C- '^i is the AWI mass with two degenerate quarks 
of bare mass mj. 

In our lattice simulations we calculate the l.h.s. of Eq. ( p!OD as a function of mi = m2 and 
m3, and extract the various constants using the following procedure. In the first step of the 
analyses we remove the contribution of the equation-of-motion operators by extrapolating 
the l.h.s to mi = 0, for fixed m^. The ratio Xq = Z^^^/Z^ Zq is then given by the intercept 
of a linear fit in m,3/2, while the slope gives Xc){bso — bo)- By choosing operators with 
different Dirac structures we are able to extract all the on-shell improvement constants, as 
well as Za, Zy and Zp/Zs- The only exception is 6t, which as discussed in Ref. ||, requires 
keeping mi 7^ m2- 

This analysis ignores 0{a^) terms. Since these can give rise to a quadratic dependence 
on quark mass, it is important to check that linear fits are adequate. In cases where the 
statistical quality of the data is good we compare linear and quadratic fits. Another check 
on the importance of 0{a?) terms is to repeat the fits using the mass ms instead of ms. In 
this case the ratio of slope to intercept gives bso — bo, which we can then compare to the 
results for bso — bo using Eq. (^). This comparison is non-trivial since the Ola"^) effects 
are different in the two cases. We stress, however, that, unless otherwise stated, the results 
presented below are from fits with respect to m.3. 

We note that, up to this point in the analysis, we do not need to introduce the off-shell 
improved operators. When we send m-i = m-2 — > we are removing the contact term between 
P and O |, and so on-shell improved operators suffice. 

This is no longer true, however, in the second step of our analysis. Here we keep m-i = m.2 
non-zero, so the contact term remains. We determine the linear combination c'p + c'q from 
the slope, so, of a linear fit of the l.h.s. of Eq. (1^) with respect to ami at fixed ams. In 
this way, for each afhs, we obtain the estimate 

c'p + c'a = 2so - Xo (bso "bo- ^bA) ■ (14) 

By choosing O = S,P,A,V,T we can determine all five c'q. Details of this part of the 
calculation are presented in Sec. pCll. 



III. SIMULATION PARAMETERS 

The parameters used in the three sets of simulations are given in Tab. |. The table also 
gives the labels used to refer to the different simulations in the following. For the lattice 



Label 


/? 


csw 


a-i (GeV) 


Volume 


L (fm) 


Confs. 


X4 


60TI 


6.0 


1.4755 


2.12 


16^ X 48 


1.5 


83 


4-18 


60NPf 
60NPb 


6.0 


1.769 


2.12 


16^ X 48 


1.5 


125 
112 


4-18 
27-44 


62NP 


6.2 


1.614 


2.91 


24^ X 64 


1.65 


70 
70 


6-25 
39-58 



TABLE I: Simulation parameters, statistics, and the time interval in x^ defining the volume V 
over which the chiral rotation is performed in the AWL The source J is placed at t = 0. 

scale a we have taken the value determined in ReL |T^ using ro, as it does not rely on the 
choice of the fermion action for a given (3. In this study what we mostly need is the change 
in scale, a(/3 = 6.2)/a(/3 = 6.0) ~ 0.73, which is much less sensitive to the physical quantity 
used to set a. 

In Tab. || we give the values of the hopping parameter k we use, along with the cor- 
responding results for am and aM^^. We also quote three estimates of Kc, obtained using 
quadratic fits, corresponding to (1) the zero of m with mass independent ca (see Sec. 0), (2) 
the zero of m with chirally extrapolated c^, and (3) the zero of M^. These are labeled Kc , 
Kc , and Kc respectively. In this paper we use Ki henceforth and drop the superscript. 

For each set of simulation parameters the quark propagators are calculated using Wup- 
pertal smearing [15|. The hopping parameter in the 3-dimensional Klein-Gordon equation 



used to generate the gauge-invariant smearing is set to 0.181, which gives mean squared 
smearing radii of {ra)'^ = 2.9 and 3.9 for f3 = 6.0 and 6.2 respectively. 

For the 60NP data set we have investigated the dependence of our results on the time 
extent of the region of chiral rotation. As shown in Tab. |, one region (forward of the source) 
is 15 timeslices long, while the other (backward of the source) is 18 slices long. Since we 
find no significant dependence on the length of the time interval, we average the two sets 
of results (assuming statistical independence). In the 62NP calculation, we also use two 
rotation regions, this time placed symmetrically about the source, in order to improve the 
signal. 

In the 60NP data set we find two exceptional configurations. Some details of the behavior 
of the pion correlator on these configurations are discussed in Appendix |^. The effect is 
most severe at the lightest quark mass, nj. We do not discard these configurations, but we 
do neglect all data with the lightest two quark masses, i.e., the hq and nj points are not 
used in the final analyses of 60NPf and 60NPb data. In the analysis of the 60TI data 
we exclude ki and Ky since the former is too heavy and the latter may have contamination 
from exceptional configurations. 



IV. RESULTS 



We begin with some general comments concerning our analysis. First, all our quoted 
results are obtained using correlation functions at zero spatial momentum. We have nu- 
merical data for non-zero momentum correlators, which lead to consistent results but with 
larger errors. Second, we use only the diagonal part of the covariance matrix when fitting 
the time dependence of correlators, or of ratios of correlators. Fits using the full covariance 
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FIG. 1: Estimates of Kc by extrapolating 62NP data for m and M^. We show quadratic fits to fh 
for the two cases discussed in text (octagons label points with mass-dependent ca and pluses label 
points with chirally extrapolated ca), and a quadratic fit to M^ (diamonds). 



matrix (which incorporates the correlations between timeslices) were not, in general, stable. 
Where we could perform such fits, we found results within la of those presented. Finally, 
fits to the quark mass dependence are also done ignoring correlations between the points at 
different masses, since our statistics are insufficient to include them. Because of the latter 
two comments, we can make no quantitative statement about goodness of fit. Nevertheless, 
assuming that the fits are good, the errors in the fit parameters, which are obtained using 
the Jackknife procedure, should be reliable. 

We begin with our results for Kc, which is needed to define the vector Ward identity 
(VWl) quark mass m. To determine k^, we make a quadratic fit of the AWl mass, rh, 
and M^ versus the tree-level quark mass parameter 1/2k. Fits to fh include only degener- 
ate quark combinations as it simplifies Eq. (P?]). Fits to M^ include both degenerate and 
non-degenerate combinations as they do not show any noticeable dependence on the mass 
difference. For the non-degenerate cases we define 2/kjj = 1/kj + 1/kj. An example of the 
resulting fits is shown in Fig. |l|. The estimate of Kc from fh should be the same whether 
we use the mass-dependent value for ca or the chirally extrapolated value in Eq. (plSf) (see 
Sec. 0). As evident from Fig. |l], the quality of both these fits is very similar and the two 
values are consistent. 

Our results for 1/kc from quadratic fits to M^ are significantly smaller than those from 
fits to m. To highlight this discrepancy we present, in Tab. |I|, the values of aM^^ at the 





60TI 


60NP 




62NP 


Label 


K 


am 


aMj, 


K 


am 


aMj^ 


K 


am 


aMyr 


Ki 


0.11900 


0.443(8) 


1.530(1) 


0.1300 


0.144(1) 


0.711(2) 


0.1310 


0.1345(6) 


0.609(1) 


^^2 


0.13524 


0.105(1) 


0.571(2) 


0.1310 


0.118(1) 


0.630(2) 


0.1321 


0.1054(4) 


0.522(1) 


K3 


0.13606 


0.084(1) 


0.504(2) 


0.1320 


0.092(1) 


0.544(2) 


0.1333 


0.0727(3) 


0.418(1) 


K4 


0.13688 


0.063(1) 


0.431(2) 


0.1326 


0.075(1) 


0.488(2) 


0.1339 


0.0560(2) 


0.360(2) 


K5 


0.13770 


0.042(1) 


0.348(3) 


0.1333 


0.056(1) 


0.416(2) 


0.1344 


0.0419(2) 


0.307(2) 


Kg 


0.13851 


0.020(1) 


0.244(4) 


0.1342 


0.032(1) 


0.308(3) 


0.1348 


0.0306(2) 


0.261(2) 


K-j 


0.13878 


0.013(1) 


0.195(8) 


0.1345 


0.025(4) 


0.262(12) 


0.1350 


0.0248(1) 


0.235(2) 


.^'^ 


0.13926(2) 





0.082(15) 


0.13532(3) 





0.083(20) 


0.135861(5) 





0.066(10) 


^> 


0.13925(2) 





0.086(15) 


0.13530(1) 





0.106(16) 


0.135862(4) 





0.073(09) 


fliQ 


0.13934(4) 






0.13541(3) 






0.13594(2) 







TABLE II: Values of the hopping parameter used in the various simulations, and the corresponding 
pseudoscalar mass aM-n- and quark mass am defined using the mass-dependent ca (see Sec. 0). 
The three estimates of Kc, obtained using quadratic fits, correspond to (1) the zero of fh with mass 
dependent ca, (2) the zero of fh with chirally extrapolated ca, and (3) the zero of M^. We quote 
the extrapolated value of uMt,^ for cases (1) and (2). 



Kc determined from fits to fh. Such a discrepancy has been observed previously (see, e.g., 
Ref. |16D, and can be attributed to a combination of quenched chiral logarithms (the effect of 
which is to cause M^ to curve downward at small quark masses [|17[ |l8l ) and chiral symmetry 



breaking by the action (which allows aMT^ifh = 0) oc a^/^ and a? effect, respectively, for 
tadpole- improved and non-perturbatively improved actions). These contributions can, in 
principle, be distinguished by the behavior of the intercept aM7r(m = 0). Quenched chiral 
logarithms are a continuum effect, implying that the intercept should be the same for 60 
TI and 60NP simulations, and that it should scale roughly proportional to a. By contrast, 
explicit chiral symmetry breaking implies a reduction in the intercept when going from 60TI 
to 60NP data sets, and an a? scaling in NP simulations. In our analyses this latter effect 
is expected to be small since fh is determined by a fit over a large range of time slices 
where the pion dominates. If these fits had extended to t = 00, then fh and M^ would 
necessarily vanish at the same k. Our results are consistent with the dependence expected 
from quenched chiral logarithms. The large residual M^, therefore, points to the need to 
include the effect of quenched chiral logarithms in the extrapolation. 

From the fit fh versus 1/2/t; we also obtain the combination {hp — 6a + &m) using Eq. (P?]). 



This is discussed in Sec. XI 



In Tabs. |Tl|, |I^, 0, and 0, we collect our results from the various Ward identities, except 



for estimates of cy which are given in Tab. [VI I| . Each identity allows us to extract one or 
more combinations of on-shell improvement and normalization constants. The details of each 
of these extractions are discussed in subsequent sections. From these results, we construct 
our best estimates for the individual constants, and these are collected in Tab. [VIII| . We 
quote both a statistical error (obtained by single elimination jackknife, in which we repeat 
the entire analysis on each jackknife sample), and an estimate of the uncertainty in the 
constants due to O(a^) errors. The latter is obtained by comparing results using values 
of Ca deduced using 2-point and 3-point discretizations of derivatives, as discussed in the 





reference 


2-point 


3-point 


CA 


Eq-dD 


-0.022(06) 


-0.023(09) 


Z'v 


Eq.(17) 


+0.747(01) 


+0.747(01) 


bv 


Eq.(O) 


+1.436(27) 


+1.455(28) 


Z'v 


Eq.(|lD 


+0.747(01) 


+0.747(01) 


bv 


Eq.(17) 


+1.534(24) 


+1.535(24) 


zyz\z^. 


Eq.(Il,D 


+1.068(13) 


+1.056(14) 


4 


Eq-dD 


+0.755(06) 


+0.759(06) 


6a -bv 


Eq.(22) 


-0.513(91) 


-0.477(95) 


Z'v 


Eq.(22) 


+0.756(06) 


+0.760(06) 


bA - bv 


Eq.(22) 


-0.488(85) 


-0.452(89) 


Z'yliZl? 


Eq.(25) 


+1.207(15) 


+1.196(16) 


bA ~ bv 


Eq.(25) 


-0.668(216) 


-0.566(222) 


zl 


Eq.@|5|) 


+0.791(07) 


+0.787(07) 


zliz\zl 


Eq.dl) 


+1.029(10) 


+1.026(13) 


z%iz\z% 


Eq.(27) 


+1.066(14) 


+1.054(15) 


bp -bs 


Eq.(|6D 


-0.070(88) 


-0.055(89) 


CT 


Eq.® 


+0.087(15) 


+0.099(18) 


bA-bp + bs/2 [cA{m)] 


Eq.(27) 


+0.739(66) 


+0.703(75) 


bA-bp + bs/2 [cAm 


Eq.® 


+0.879(64) 


+0.021(46) 


bp — bA 


Eq.(31) 


-0.126(58) 


-0.125(81) 


bs -by - 2{bp - bA) 


Eq.® 


-0.588(274) 


-0.266(380) 



TABLE III: Results for the 60TI data set. 



following section. Another estimate of O(a^) errors is obtained by comparing our results 
to the previous estimates of the ALPHA collaboration 0, ^ |^ which we also include in 
Tab. VIII| along with the one-loop perturbative results discussed in Appendix 0. We quote 
both by, bA and by, bA to simplify comparison with previous results. 

We collect separately, in Tab. |^, our results for the improvement constants c'x, the 
coefficients of the equation-of-motion operators. These are discussed in Sec. [XI1| . 

The assiduous reader will notice that our results for the 60TI data differ slightly from 
those presented in Ref. §. This is for two reasons. First, we use a new method for determining 
cy. This leads to a much more precise result, and affects several other analyses which are 
dependent on cy. Second, we have made several minor improvements in our analysis, e.g. 
using quadratic instead of linear fits versus quark mass where appropriate. The set of 
configurations has not changed. 

We now discuss the salient features of our final results from Tab. |V114 Perhaps the most 
important issue is the comparison with the results by the ALPHA collaboration. Because 
we use different improvement conditions, the results for the Zx can differ by ~ a^Ag^^, 
while those for the cx and bx can differ by ~ aKqco- Numerically these are about 0.02 and 
0.15, respectively, at /? = 6.0, and 0.01 and 0.1, respectively, at /5 = 6.2. There are some 
quantities, however, where these differences can be enhanced. For example, in correlators 





reference 


2-point 


3-point 


CA 


Eq-dD 


-0.037(04) 


-0.045(07) 


Z'v 


Eq.(17) 


+0.770(01) 


+0.769(01) 


bv 


Eq.(|T7l) 


+1.429(20) 


+1.466(24) 


Z'v 


Eq.(|l7D 


+0.769(01) 


+0.768(01) 


bv 


Eq.(17) 


+1.524(14) 


+1.525(14) 


zyz\z^. 


Eq.(l4l9D 


+1.067(09) 


+1.041(13) 


^^y^ 


Eq.(^ 


+0.773(04) 


+0.775(04) 


bA -bv 


Eq.(22) 


-0.231(47) 


-0.179(57) 


Z'v 


Eq.(22) 


+0.774(03) 


+0.776(04) 


bA - bv 


Eq.(^ 


-0.216(43) 


-0.165(53) 


zVizlY 


Eq.(25) 


+1.197(09) 


+1.185(10) 


bA ~ bv 


Eq.(25) 


-0.193(91) 


-0.180(107) 


zl 


Eq.(pq,p5|) 


+0.808(03) 


+0.800(03) 


zliz\zl 


Eq.(^ 


+1.048(09) 


+1.035(11) 


zyz\z% 


Eq.(27) 


+1.049(08) 


+1.026(11) 


bp -bs 


Eq.(^ 


-0.013(55) 


+0.019(57) 


CT 


Eq.(|30D 


+0.063(07) 


+0.092(11) 


bA-bp + bs/2 [cA{m)] 


Eq.(27) 


+0.609(31) 


+0.570(53) 


bA-bp + bs/2 [cAm 


Eq.® 


+0.883(32) 


-0.052(22) 


bp — bA 


Eq.(^ 


-0.079(54) 


-0.031(74) 


bs -by - 2{bp - bA) 


Eq.(^ 


-0.331(201) 


+0.112(338) 



TABLE IV: Results for the 60NPf data set. 



dominated by the pion, contributions proportional to aBT^ = aM^/2rh, while formally of 
©(oAqcd), can be numerically much larger. These cases are discussed in more detail in the 
following sections. 

Given these estimates of the uncertainties, we find that, at /? = 6.2, there is complete 
consistency between our results and those from the ALPHA collaboration. Indeed, the 
only statistically significant difference is for Zy, which is calculated very precisely, but this 
difference is consistent with being an ~ o^Aq^^ effect. 

Moving to /3 = 6, we see that there are statistically significant differences not only for 
but also for ca and cv- For Zy the differences are consistent with the estimates of 
discretization errors given above. The difference for cy (ca) is about two (three) times the 
expected size of ~ 0.15 — this could be an enhanced 0{a) correction or an effect of higher 
order in a. Either way, what is clear is that, within 0{a) improvement, non-perturbative 
estimates of the cx have substantial uncertainties at /5 = 6. The only definite conclusion 
that we can draw is that the cx, which are zero at tree level, are small. 

We find that the various constants show a strong dependence on the value of csw- The 
relatively small change from the non-perturbative value at /3 = 6 to the tadpole-improved 
value leads to noticeable changes in most of the constants. 

One of the most surprising results of Ref. || was the large magnitude of 6^ — 6^ ~ 0.5 at 
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reference 


2-point 


3-point 


CA 


Eq-dD 


-0.036(05) 


-0.043(08) 


Z'v 


Eq.(17) 


+0.770(01) 


+0.769(01) 


bv 


Eq.(O) 


+1.424(17) 


+1.464(23) 


Z'v 


Eq.(|lD 


+0.769(01) 


+0.768(01) 


bv 


Eq.(17) 


+1.522(11) 


+1.523(11) 


zyz\z^. 


Eq.(Il,D 


+1.068(10) 


+1.041(14) 


4 


Eq-dD 


+0.766(04) 


+0.766(04) 


6a -bv 


Eq.(22) 


-0.288(43) 


-0.256(53) 


Z'v 


Eq.(22) 


+0.768(04) 


+0.768(04) 


bA - bv 


Eq.(22) 


-0.267(40) 


-0.236(50) 


Z'yliZl? 


Eq.(25) 


+1.204(11) 


+1.194(13) 


bA ~ bv 


Eq.(25) 


+0.007(106) 


+0.043(126) 


zl 


Eq.@|5|) 


+0.806(04) 


+0.797(04) 


zliz\zl 


Eq.dl) 


+1.061(10) 


+1.050(14) 


z%iz\z% 


Eq.(27) 


+1.051(08) 


+1.027(12) 


bp -bs 


Eq.(|6D 


-0.114(44) 


-0.097(44) 


CT 


Eq.® 


+0.057(10) 


+0.084(13) 


bA-bp + bs/2 [cA{m)] 


Eq.(27) 


+0.596(33) 


+0.547(58) 


bA-bp + bs/2 [cAm 


Eq.® 


+0.881(33) 


-0.051(23) 


bp — bA 


Eq.(31) 


-0.058(54) 


+0.002(81) 


bs -by - 2{bp - bA) 


Eq.® 


-0.379(247) 


+0.096(439) 



TABLE V: Results for the 60NPb data set. 



/3 = 6 with tadpole- improved csw- This difference is predicted to be very small (0.002) in 



1-loop perturbation theory, and even assuming the 2-loop term to be 



a: 



suggests a much 



smaller value ~ 0.02 (0.015) at /3 = 6 (6.2). We find that the measured difference is reduced 
to ~ 0.3 using the non-perturbative csw, and further reduced to ~ 0.1 at /5 = 6.2. While 
the latter difference is small enough to be accounted for by the expected qAqcd uncertainty, 
the larger result at /5 = 6 may indicate higher order uncertainties. 

The other differences between 6's are more stable, and are consistent with perturbative 
predictions within the 0{a) uncertainties. The same is true of our final results for the 6's 
themselves; the largest difference is for by and is ~ 2aAQCD- 

In fact, allowing for (1 — 2)aAQCD discretization uncertainties, the only non-perturbative 
result which is in disagreement with perturbation theory is Zp/Zg. A very large two loop 
effect, ~ — 4q;^, is required to bring the results into agreement. This finding is consistent 
with those of the APE collaboration who argue that Zp — 1 is significantly underestimated 
by 1-loop perturbation theory |TP . 



Concerning the statistical errors, we see a substantial improvement in the signal between 
(3 = 6.0 and 6.2. It is also noteworthy that the errors in our estimates are comparable to 
those from the ALPHA collaboration. While a precise comparison of efficacies is difficult 
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reference 


2-point 


3-point 


CA 


Eq-dD 


-0.032(03) 


-0.038(04) 


Z'v 


Eq.(17) 


+0.787(00) 


+0.787(00) 


bv 


Eq.(O) 


+1.304(10) 


+1.312(10) 


Z'v 


Eq.(|lD 


+0.787(00) 


+0.787(00) 


bv 


Eq.(17) 


+1.422(08) 


+1.422(08) 


zyz\z^. 


Eq.(Il,D 


+1.091(05) 


+1.084(05) 


4 


Eq-dD 


+0.788(02) 


+0.790(02) 


6a -bv 


Eq.(22) 


-0.111(27) 


-0.071(28) 


Z'v 


Eq.(22) 


+0.788(02) 


+0.791(02) 


bA - bv 


Eq.(22) 


-0.109(26) 


-0.071(27) 


Z'yliZl? 


Eq.(25) 


+1.185(04) 


+1.181(05) 


bA ~ bv 


Eq.(25) 


-0.092(62) 


-0.115(59) 


zl 


Eq.@|5|) 


+0.818(02) 


+0.813(02) 


zliz\zl 


Eq.dl) 


+1.085(04) 


+1.077(05) 


z%iz\z% 


Eq.(27) 


+1.077(05) 


+1.071(05) 


bp -bs 


Eq.® 


-0.086(23) 


-0.075(23) 


CT 


Eq.® 


+0.051(07) 


+0.078(07) 


bA-bp + bs/2 [cA{m)] 


Eq.(27) 


+0.626(24) 


+0.619(29) 


bA-bp + bs/2 [cAm 


Eq.® 


+0.850(19) 


+0.123(17) 


bp — bA 


Eq.(31) 


-0.086(26) 


-0.062(34) 


bs -by - 2{bp - bA) 


Eq.® 


+0.047(106) 


+0.176(137) 



TABLE VI: Results for the 62NP data set. 



because of different systematic errors, and different ensemble and lattice sizes, we conclude 
that our method is competitive. 

V. CALCULATION OF ca 

The determination of ca is central to the extraction of all quantities that are obtained 
using the axial Ward identity Eq. (|iy) since ca enters in 6S [see Eq. (H)]. Its evaluation 
uses the AWI with no operator present in the domain of chiral rotation. In particular, ca is 
adjusted so that the ratio 



EAMA, + acAd,P](^^\x, t)j(^nO)) 

E,(P(^^-)(x,t)Jo-^)(o)) 



(15) 



which defines the quark mass rhij, is independent of the source J and the time t at which it 
is evaluated. Since this criterion is automatically satisfied when the correlators are saturated 
by a single state, the determination of ca relies on the behavior of excited state contributions 
at small t. 
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62NP 




2pt 


3pt 


CA{m) 


ca{0) 


CA{m) 


ca{0) 


extrap. 


-0.115(63) 


-0.087(62) 


-0.032(64) 


-0.096(61) 


1/mfit 


-0.086(15) 


-0.102(17) 


-0.172(20) 


-0.123(19) 


slope ratio 


-0.094(19) 


-0.094(19) 


-0.107(19) 


-0.109(19) 


60NPf 




2pt 


3pt 


CA{m) 


ca(0) 


CA{m) 


ca{0) 


extrap. 


-0.094(56) 


-0.060(57) 


+0.046(68) 


-0.048(63) 


1/mfit 


-0.131(26) 


-0.205(38) 


-0.363(71) 


-0.209(53) 


slope ratio 


-0.116(19) 


-0.116(19) 


-0.113(26) 


-0.120(26) 


60NPb 




2pt 


3pt 


CA(m) 


ca(0) 


CA{m) 


ca{0) 


extrap. 


-0.119(78) 


-0.086(78) 


+0.013(87) 


-0.067(81) 


1/mfit 


-0.071(38) 


-0.157(45) 


-0.359(86) 


-0.171(63) 


slope ratio 


-0.097(30) 


-0.096(31) 


-0.092(37) 


-0.102(37) 


60TI 




2pt 


3pt 


CA{m) 


ca{0) 


CA{m) 


ca{0) 


extrap. 


-0.483(124) 


-0.468(125) 


-0.337(135) 


-0.451(131) 


1/mfit 


-0.143(63) 


-0.162(65) 


-0.367(80) 


-0.158(71) 


slope ratio 


-0.253(48) 


-0.252(49) 


-0.222(53) 


-0.244(55) 



TABLE VII: Results for cy- See text (sec. VII) for details 



To implement Eq. ([TsD one has to choose how to discretize the derivatives. Note that 
all choices lead to the same improvement and normalization constants at the order we are 
working, i.e. up to 0{a) and O(a^) errors, respectively. This is because the difference 
between discretizations is explicitly proportional to a^. Thus investigating the sensitivity to 
the choice of discretization gives information on the size of higher order discretization errors. 

We limit our study of this issue to the comparison between two discretization schemes. 
Both are based on a mixture of 2-point and 3-point discretizations. This terminology is 
explained in Ref. |, and is exemplified by dxf{x + 0.5) — > (/(a; + 1) — f{x))/a (2-point) 
and dxf{x) -^ {f{x + 1) — f{x — l))/2a (3-point). Results from both schemes are quoted in 
Tabs. 0,0 and |VI[ 

In our first scheme, we implement Eq. (p!5| ) using 2-point discretization. In the subsequent 
calculations, based on the AWI of Eq. (|T^) we use the same 2-point discretization in 6S 
[Eq. (^] as in Eq. (P^, and replace the continuum integral by a simple sum. For the 
derivatives within the operators O and 60, however, we use 3-point discretization. In our 
second scheme, we repeat the calculations using the value of Ca obtained when enforcing 
Eq. (|15D with a 3-point discretization for the derivatives. The remainder of the calculation 
is done with the same discretizations for 6S, O and 50 as in the 2-point scheme but the 
new value for ca- 

There is a subtlety in the comparison between results from the two schemes. It follows 



13 





/3 = 6.0 


13 = 6.2 




LANL 


LANL 


ALPHA 


P. Th. 


LANL 


ALPHA 


P. Th. 


csw 


1.4755 


1.769 


1.769 


1.521 


1.614 


1.614 


1.481 


A 


+0.747(1) 


+0.770(1) 


0.7809(6) 


+0.810 


+0.7874(4) 


+0.7922(4) (9) 


+0.821 


z\ 


+0.791(7)(4) 


+0.807(2)(8) 


0.7906(94) 


+0.829 


+0.818(2)(5) 


+0.807(8)(2) 


+0.839 


zyzi 


+0.811(9)(5) 


+0.842(5)(1) 


N.A. 


+0.956 


+0.884(3)(1) 


N.A. 


+0.959 


CA 


-0.022(6)(1) 


-0.037(4)(8) 


-0.083(5) 


-0.013 


-0.032(3)(6) 


-0.038(4) 


-0.012 


cy 


-0.25(5)(3) 


-0.107(17)(4) 


-0.32(7) 


-0.028 


-0.09(2)(1) 


-0.21(7) 


-0.026 


CT 


+0.09(2)(1) 


+0.06(1)(3) 


N.A. 


+0.020 


+0.051(7)(17) 


N.A. 


+0.019 


~bv 


+1.44(3)(2) 


+1.43(1)(4) 


N.A. 


+1.106 


+1.30(1)(1) 


N.A. 


+1.099 


by 


+1.53(2) 


+1.52(1) 


+1.54(2) 


+1.273 


+1.42(1) 


+1.41(2) 


+1.254 


bA - by 


-0.51(9)(4) 


-0.26(3)(4) 


N.A. 


-0.002 


-0.11(3)(4) 


N.A. 


-0.002 


bA -by 


-0.49(9)(4) 


-0.24(3)(4) 


N.A. 


-0.002 


-0.11(3)(4) 


N.A. 


-0.002 


bp -bs 


-0.07(9)(2) 


-0.06(4)(3) 


N.A. 


-0.066 


-0.09(2)(1) 


N.A. 


-0.062 


bp -bA 


-0.126(58)(1) 


-0.07(4)(5) 


N.A. 


+0.002 


-0.09(3)(3) 


N.A. 


+0.002 


bA 


+0.92(10)(6) 


+1.17(4)(8) 


N.A. 


+1.104 


+1.19(3)(5) 


N.A. 


+1.097 


bA 


+1.05(9)(4) 


+1.28(3)(4) 


N.A. 


+1.271 


+1.32(3)(4) 


N.A. 


+1.253 


bp 


+0.80(11)(6) 


+1.10(5)(13) 


N.A. 


+1.105 


+1.11(4)(7) 


N.A. 


+1.099 


bs 


+0.87(14)(4) 


+1.16(6)(11) 


N.A. 


+1.172 


+1.19(4)(6) 


N.A. 


+1.161 



TABLE VHL Final results for improvement and renormalization constants. The first error is 
statistical, and the second, where present, corresponds to the difference between using 2-point and 
3-point discretization of the derivative used in the extraction of ca- 



from the relation 

(9^[(4'-^°^"*))^(c^) - 2mP(2-Pomt)](^ + a/2) J(0)) + 

(9^[(Af-^""*))^(c^) - 2mP(2-Pom<)](^ _ a/2)J(0)) = 

2(a^[(Af-^°^"*))^(cA - am/2) - 2mP](t) J(0)) + 0{a') . (16) 

that the Ola"^) differences between 2-point and 3-point discretizations can be absorbed by 
shifting Ca —^ Ca — am/2 in the latter scheme. Thus, if one were to fit to the same range of 
timeslices with appropriate weights, as defined by Eq. (|TB|), the difference between ca from 
2-point and 3-point determinations would be of 0{a^) in the chiral limit. This difference 
would then not be useful as an indicator of 0{a) discretization errors. 

In practice, however, our fits do not weight the points appropriately for the relation (p!6D 
to be relevant. In particular, we find that using the 2-point scheme, the best fits are for 
t > 2 relative to the source at t = 0, where t = 2 (which corresponds to evaluating the 
derivative at t = 2.5) is the earliest timeslice at which there are no contact terms for either 
discretization scheme. On the other hand, for the 3-point scheme, we are not able to include 
the point at t = 2 as the O^a"^) errors are too large and the fit has poor quality (this was 
checked by turning on the full covariance matrix). Because of this, the resulting values of 
Ca do differ at 0{a), and we take this the difference as an estimate of the size of the higher 
order discretization errors. 



14 





60TI 


60NPf 


60NPb 


62NP 


C'y + dp 

c'a + c'p 
2dp 

c's + c'p 

drp + c'p 


+2.75(23) 

+2.30(46) 

-1.96(152) 

+2.02(21) 

+2.26(33) 


+2.82(15) 
+2.43(24) 
+0.88(97) 
+2.44(13) 
+2.40(18) 


+2.68(19) 
+2.12(31) 
-0.65(57) 
+2.40(13) 
+2.27(20) 


+2.62(8) 
+2.43(14) 
+1.82(24) 
+2.40(7) 
+2.42(9) 


4 

c'p 
c's 

C-rp 


+3.72(73) 
+3.28(94) 
-0.98(76) 
+3.00(73) 

+3.24(75) 


+2.38(50) 
+1.99(56) 
+0.44(49) 
+2.00(48) 
+1.96(49) 


+3.00(37) 
+2.45(46) 
-0.33(29) 
+2.72(33) 
+2.60(38) 


+1.72(16) 
+1.53(20) 
+0.91(12) 
+1.49(14) 
+1.51(15) 



TABLE IX: Results for ofF-shell mixing coefficients. 



In our final compilation, Tab. |V11I| , the central values are from the 2-point discretization, 
while the difference between the two discretizations is quoted as a systematic error. We note 
that the ALPHA collaboration has used 3-point discretization of all derivatives. This does 
not, however, imply that their results should be more closely comparable to ours based on 
the ca with 3-point discretization, since there are other differences in the calculations. 

To use Eq. ( |T5|) we must also choose the source J. Different sources produce different 
admixtures of the ground and excited states, and thus have varying sensitivities for deter- 
mining Ca- Furthermore, different sources give values for ca differing by 0{a) (or 0(1) if the 
action is not fully 0{a) improved). We have investigated source dependence using results 
from a separate calculation performed on 170 quenched lattices of size 32^ x 64 at /3 = 6.0 
using the Wilson {csw = 0) 12^] and tadpole- improved clover {csw = 1-4785) |^ actions. 
(The slightly different value of csw = 1-4755 used in the 60TI calculation was an oversight.) 
The results from three different sources are shown in Fig. ^. The sources are J = A^ and 
J = P, both with wall source smearing, and J = P with Wuppertal smearing. We do not 
present the J = A^ data with Wuppertal smearing as that correlator is dominated by the 
ground state already at t ~ 4, and is thus very insensitive to ca- Results from the Wilson 
action depend substantially on the source, even in the chiral limit. This is as expected since 
the action is not improved, leading to large, 0(1), variations in ca- Also, as expected, there 
is a marked convergence when using the tadpole- improved action. Indeed, results from the 
three sources are consistent within errors (and linear extrapolation to the chiral limit gives 
a result, ca = —0.026(2), consistent with our 60TI result quoted in Tab. |m| ), and have 
similar sensitivity in determining ca- Because of this, we have chosen to use only J = P 
with Wuppertal smearing in the simulations devoted to calculating improvement constants- 

We illustrate our determination of ca (with 2-point discretization) using the non- 
perturbatively improved action in Figs- ^ and |. We tune ca so as to extend the plateau 
to the earliest timeslice t = 2 at which there are no contact contributions (the source is 
at t = 0). We have enough sensitivity to clearly distinguish ca from zero- At /5 = 6 we 
can also distinguish ca from that obtained by the ALPHA collaboration for the chiral limit 
{ca = —0-083(5))- This difference remains after we extrapolate our results to the chiral 



15 



0.05 



o 



-0.05 



T 1 1 r 



1 1 1 r 



~i 1 1 r 



w <> 



<> 



<> 



<> 



X J=A4, Wall 
o J=P, Wupp 
o J=P, Wall 



<> 



<> 








0.05 0.1 



0.15 



FIG. 2: Estimates of ca versus the quark mass for three different sources J as discussed in the text. 
For the Wilson action (W), estimates of ca from the three J are very different. The improvement 
in going to the tadpole-improved clover action (TI) is dramatic, and the three sets of data collapse 
together. We show a linear fit to this combined tadpole-improved clover data. 

limit (giving ca = — 0.037(4)^for the two-point discretization and ca = 0.045(7) for the 
three-point discretization). At /? = 6.2 our results for ca differ from the ALPHA value, 
Ca = —0.038(4), at non-zero quark mass (as shown in the Figs. § and ^, but after chiral 
extrapolation they are consistent with the ALPHA result. This extrapolation (which is done 
using a linear fit to the masses K2 — K5) is shown in Fig. |[ 

In our previous paper [§] we used tadpole- improved fermions at /3 = 6, and found a result 
inconsistent with that of the ALPHA collaboration, as can be seen from Tab. IVIIIL We did 



not, however, have enough information to determine the source of this difference. Our new 
result shows that while increasing csw to its non-perturbative value moves ca towards the 
ALPHA result, a significant difference of ~ 0.046 remains. This difference is presumably due 
to higher order discretization errors. It is striking, however, that the difference is reduced 
substantially by changing (3 from 6.0 to 6.2. Since a^ only halves between /3 = 6 and 6.2, 
this suggests that even higher order discretization errors are playing a dominant role. By 
contrast, the reduction in the difference between our results for 2- and 3-point discretizations 



^ The difference in ca in the 60NPf and 60NPb estimates is due only to the different number of configu- 
rations analyzed. Since the 60NPb sample is a subset of 60NPf, we quote the 60NPf result as our best 
estimate 



16 



0.2 



C\2 



0.15 



0.1 



X 



X 



X 



(D 



. ^MMXM . y ... ii!^ .,, X .... >if ... X ... X"^ ' "'y ' "^ ' "y"X .... ^ ... )i; .. .X ... X .... y ,, 

4> ... ^„ , ... 4) ... o .., |^ .., ^ .... i^ .. „^j^„„^ ... ^ ... ^ ... a? .... <j) .... ^ . „ 



<1> <I> 



<5> 



<3> 



-O 



J I I L 



X Ca= 

o Ca=-0.022 

o Ca=-0.083 



J I I L 







10 



20 



FIG. 3: Estimates of 2mij for different values of ca illustrated using i = j = K3 in the 60NPf 
data set and 2-point discretization. For this value of quark mass, setting ca = —0.022 extends the 
plateau to the earliest time slice t = 2 at which there are no contact contributions. The fit for 
CA = —0.083, the value obtained by the ALPHA collaboration in the chiral limit, and ca = are 
included to illustrate sensitivity. 



is consistent with being an a^ effect. 

It is interesting to compare our non-perturbative results for ca with perturbative es- 
timates. We see from Tab. |VIII| that the 1-loop result (~ 0.2 x a) gives a substantial 



underestimate. To explain the difference one needs a large two-loop term, ~ a , which, 
using the values quoted in Appendix ^, is 0.018 and 0.016 for (3 = 6 and 6.2, respectively. 

We close with a comment on the practical implementation of the AWL To the accuracy 
we are working, we can use, in SS, either the appropriate mass- dependent ca or its value 
in the chiral limit. We prefer the former, and use it throughout, because it maintains the 
relation df^{Aj)fi — 2m.ijP^'''^^ = at finite quark masses on the states used to tune ca- Our 
data suggest that this relation receives only small corrections on other states relevant to the 
AWL This ensures (for the ca obtained using 2-point derivatives) that the ratio in Eq. (p!0[) 
is nearly independent of the time slice of the insertion of the improved operator and the 
volume V of chiral rotation. We stress, however, that when the axial current appears as 
an operator in the AWI, we use the chirally extrapolated ca to give our central values (see 
Section |ID, and use the mass-dependent ca to give an indication of the size of higher order 
discretization errors. 
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FIG. 4: Estimates of 2'mij for different values of ca illustrated using i = j = k^ and the 62NP 
data set. For this quark mass, ca = —0.0209 extends the plateau to the earliest allowed time slice 
t = 2. To show sensitivity to the tuning we contrast this best fit with those using ca = and 
CA = -0.040. 



VI. Z^ AND bv 

The matrix elements of the vector charge J d^xV^ {^)^ with m2 = rris, are fixed by the 
charge of the states, and allow a determination of Zy as a function of the quark mass. Our 
best signal is for the matrix element between pseudoscalar mesons: 



Z^{1 + hvam2) 



(E..Pa2)(f,r)J(2i)(0)) 



(17) 



with r > t > and J = P or A4. The two sources have comparable signal, and the 
final results are obtained by averaging the two estimates when constructing the jackknife 
ensemble. Note that the 0{a) improvement term in Vj does not contribute. Zy and by are 
then extracted by fitting the data as a function of m2. 

As an illustration we describe the procedure for the 62NP data set. The quality of the 
data is very good, as shown in Fig. |^. A linear fit is clearly inadequate, so we use a quadratic 
fit 



Zy = 0.7874(4) [1 + I.304(10)m2a + I.062(52)(m2a)^] . 



(18) 
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FIG. 5: The chiral extrapolation of ca for 62NP data. Diamonds label all mass combinations and 
stars highlight the ten combinations of K2 — k,5 used in the fit. 



The intercept is our result for Zy, while the coefficient of the linear term, i.e. the slope in 
the chiral limit, is our result for by. Note that if we had simply used a linear fit over our 
mass range, the result for by would have been 1.469(9), in complete disagreement with our 
quoted result. 

We can also fit Zy as a function of m = 1/2k, — 1/2k,c- This provides a consistency check 



for Zy, and a direct determination of by. The fit gives 

Zy = 0.7871(3) [1 + 1.422(8)ma + 0.05(4)(ma)2] 



(19) 



In this case the quadratic term is small. The intercept is consistent, at the 2-cr level, with 
that from Eq. (|TB|). We can use these two fits to also extract the combination 



Z'p/Z'aZI) 



from the ratio of the coefficients of the linear term as explained in Sec. |X| The results are 



given in Tabs. pIHVI|, and are consistent with those obtained from the axial Ward identity, 
Eq. (p6D, even though the Ola"^) errors could have been different in the two methods. 



Our results for Zy and by are compared with those from the ALPHA collaboration in 



Tab. [Vlll] . There are small differences for Zy, 0.011(1) and 0.005(1), respectively, at /3 = 6.0 
and 6.2. These are of the expected magnitude for 0{a^) differences, and are consistent with 
0{a^) scaling. The results for by are, on the other hand, already consistent. 

The difference between 1-loop tadpole-improved perturbation theory and our non- 
perturbative Zy is 0.040(1) at /5 = 6.0 and 0.034(1) at /3 = 6.2, where only statistical 
errors have been considered. Recall that the discretization errors are expected to be of 
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FIG. 6: Linear and quadratic fit to Zy versus m2 for the 62NP data set. 



size (aAqco)^ ~ 0.02 and 0.01, respectively, while the missing two loop perturbative terms 
should be ~ a^ ~ 0.02 and 0.016, respectively. Thus the deviation from perturbation theory 
is of the expected size, and the scaling behavior is closer to O(a^) than to 0{a^). The 
numerical values are consistent with p^ 2a^. 

The non-perturbative results for by exceed the 1-loop estimates by 0.24(2) and 0.16(2), 
respectively, at the two couplings. These differences are much larger than the missing two- 
loop contributions, but are consistent with a discretization error of size ~ l.SaAqcD- 

Results for Zy are needed to calculate the vector decay constants and semi-leptonic form 
factors of D and B mesons. Note that, at /3 = 6.2, the charm and bottom quark masses 
are, in lattice units, roughly 0.5 and 2.0, respectively, to be compared to our largest mass of 
0.13. It is thus important to ascertain to what mass the fits, given in Eqs. (|l^) and (19), can 
be used reliably. To address this issue we show in Fig. ^ how the two fits extend to higher 
quark masses for (3 = 6.2. A plot of the quantity m/m [Eq. (^)] which we use to convert fha 
to ma is also included. We also show the recent non-perturbative results for Zy obtained 
by the UKQCD collaboration |2^. Comparing our fits with the UKQCD data we find that 
both fits provide reliable estimates (to within 2%) up to the charm quark mass, with the fit 
to Eq. (18) being slightly better. In fact, over the range < ma ^ 0.5, truncating Eq. (19) 
at the linear order fits the UKQCD data to within 1% as already noted by them. Beyond 
m,a ~ 0.5 the two fits start to deviate, and their validity near the bottom quark mass needs 
to be examined. 
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FIG. 7: Predictions for Zy at /? = 6.2 obtained by extending our fits, eqs. (18) and (19), to larger 
quark masses. The result for mjm is also shown, as are data points from the UKQCD cohaboration. 

VII. cyAND6^-V 



We now turn to the analyses of the various 3-point axial Ward identities, and first consider 
the determination of the improvement coefficient cy. A precise determination of cy is impor- 
tant both for phenomenological apphcations and because the uncertainty in cy contributes 
significantly to errors in Z\, Z^/Z^, ct, and c^. We have investigated several methods, and 
obtain the best results by enforcing 






dfj^i^\ 






,(20) 



where the dependence on cy enters only on the r.h.s.. We emphasize two important features 
of this method. First, it does not require knowledge of the normalization constants Za and 
Zy, since these appear in the same combination on both sides of Eq. (|20D . Second, the 
relation holds for any value of the quark masses, since the contact terms are the same on 
both sides [see Eq. (0)]. The determination of cy does, however, require knowledge of ca, 
which enters both in SS and in {Ajy^ on the l.h.s.. 
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The two correlators on the Lh.s. are dominated by the pion channel and the signal is 
excellent in the individual correlators as well as in the ratio. The latter is illustrated in 
Fig. p. On the other hand, the correlators on the r.h.s. are dominated by the ai intermediate 



state, for which the signal is not as good. We illustrate this by showing, in Figs. |^ and |T0 
the terms independent of and proportional to cy. It turns out that the difference between 
the l.h.s. and the cy independent term on the r.h.s. is about 2% of the individual terms, 
and is comparable to the error, which is dominated by that from the term on the r.h.s.. 
Nevertheless, as explained below, we can extract cy with reasonable precision. To do this it 
is convenient to rewrite Eq. (^) in terms of the following two quantities: 



T.y{5sf^ {Vi)T\y.y^) J'^'^Ko)) E,-(^^r vF'\y,y^) a?'\o)) 

T.yU,)T\y,y,) J(3i)(0)) T.y{{A,)f\y,y,) Af ^(O)) 

Y. y{5Sf^ ad,Til'\y,y,) Af\ Q)) 



N 



— " ' ^(l3)/-,„, N .(31)^"- ■ ^ ' 



such that Cy = N/ D. 

The data exhibit three interesting features: 

• Both A^ and D are, to a good approximation, linear in rhi — m^, as illustrated in 
Fig. |ll] (A^ shows a weak dependence on m^ + mi as well). 

• Close to rhi = m^ both A^ and D vanish. However, since the discretization errors in 
N and D are different, they vanish at slightly different points. As a result the ratio 
N/D is very poorly determined when fhi = rh^, and cy = N/D shows a spurious 
l/{fhi — ms) singularity, as illustrated in Fig. |l^. 

• Estimates of cy for the combination {rhi, rhj} are highly anti-correlated with those for 
{rhj,rhi}. Estimates of cy for rhi < ^3 are consistently more negative as shown in 
Fig.0 

Because of the spurious singularity mentioned above, we explore the following three 
approaches to determine cy- 

• Linearly extrapolate each of the three ratios of correlators to fhi = fri2 = 0, working at 
fixed non-zero rhs so as to avoid the singularity, and then solve for cy. The weighted 
average over the different ma points is quoted in the first row in Tab. |V11| . This method 
yields estimates with the largest uncertainty, as illustrated in Fig. |TB|. 

• Fit N/D to the form Cy + Cy /{rhi ~ ''^3) (as illustrated in Fig. 0) and use cy = Cy . 
We find that the result is insensitive to the range of quark masses used; the results 
quoted in Tab. |V11| are based on fits to k,2 — 1^5 for 60TI and 60NP and ni — kq for 
62NP. 

• Fit A^ and D separately to the form a + 7(mi — "n^a), and take cy to be the ratio of 
the slopes, ■Jn/id- This is legitimate since cy is given, in principle, by N/D for all 
quark masses. This method avoids the use of the intercepts, a^ and ao, which, being 
small, have larger discretization errors. 
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FIG. 8: Illustration of the quality of the signal for the l.h.s. of Eq. (|2 
all four cases all quark propagators correspond to H3. 



for the four data sets. In 



For each of these methods we evaluate cy for four variants of c^: for both the usual choices 
of 2-point versus 3-point discretization of (94^4 when determining ca using Eq. (0), we use 
mass dependent and chirally extrapolated values of ca in the operator {Ajy^ appearing in 
the denominator on the l.h.s. of Eq. (pO]).^ Results are quoted in Tab. |V11| . We find that 
only for the "slope-ratio" method do all four choices for ca lead to consistent results. We 
also note that the estimates using all three methods are consistent if we use the 2-point ca 
but not for the 3-point ca- Thus we take for our best estimate the value obtained with the 
"slope-ratio" method and the 2-point (chirally extrapolated) ca- 

Our final results are collected in Tab. |VII1| . Our main conclusion is that cy, which is 
zero at tree level, remains small in magnitude. We note that although our non-perturbative 
estimate is smaller than those of the ALPHA collaboration, the difference is consistent with 
being due to aAqcd corrections. 

We have tried several other methods for determining cy- One can demand that the r.h.s. 
of Eq. (|20|) be independent of 2/4. This turns out to be roughly true for the individual ratios, 
and thus holds independent of cy- We have also tried different sources, e.g. O = Ai, J = Vi 
at zero momentum for the l.h.s. and O = A4, J = V4 at non-zero momentum for the r.h.s.. 



As noted in Sec. [y|, we always use the mass dependent ca in SS. 
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FIG. 9: Illustration of the quality of the signal for the first term on the r.h.s. of Eq. ( pO[ ) for the 
four data sets. In all four cases all the quark propagators correspond to k^. 



making the intermediate state a vector meson. We have also implemented the method of 
the ALPHA collaboration ^ in which the l.h.s. with O = Ai, J = ^ at zero momentum is 
equated to unity in the chiral limit, making use of previously determined results for Z^ /Zy 
and Zy. In all the cases we have considered, however, the final estimates have larger errors 
than those quoted above. It is noteworthy, and perhaps surprising, that our best method 
involves an intermediate axial-vector state, rather than a vector meson. 

The errors in our final result for cy are substantially smaller than those of Ref. ^ It is 
likely that part of the explanation for this improvement is our use of a different AWI and 
fitting method. 



To extract 6a — by we use the l.h.s. of Eq. (pO]), so as to avoid dependence on cy, and 
follow the procedure outlined in Sec. |T[ After extrapolating to fhi = m2 = 0, the ratio 
should be described by 



^0(1 + bAams/2) 



Z\ 



Z^{l + hvam^/2) 



(22) 



The slope with respect to m^/2 {m^/2) gives our best estimate for h^ — by {pA — by) and 



the intercept gives a second estimate of Zy. As shown in Tabs. p^[V^ , the results for Zy 
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FIG. 10: Illustration of the quality of the signal for the ratio multiplying cy in the r.h.s. of Eq. (20) 
for the four data sets, using K3 propagators in all cases. 



are consistent with those from the VWI, but with somewhat larger errors. As an example 
of the fits, for the 62NP data set we find 



1 + (&yi - bv)am3/2 
1 + (^A - bv)am3/2 



1.269(3) (1 -0.111(27)am3/2) 
1.268(3) (1 -0.109(26)am3/2) 



(23) 



The quality of the fits is shown in Fig. |T^. Even though the intercept and the slope are almost 
identical, they are consistent with the expected relation {b^ — by) = {Z\Z^/Z%){bA — by) ~ 
0.92(6^ — by) within the errors. 

Since the correlators on the l.h.s. of Eq. (^) involve pion intermediate states, higher 
order discretization errors can be enhanced as noted in Sec. |V| For example, a change in 
the value of ca used in the denominator, Ac^ ~ qAqcd, leads to a change in Zy of size 
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FIG. 11: 62NP data for N and D used to extract cy and defined in the text, plotted as a function 
of mi — fh^. 

The ratio B^, = M^/fh is much larger than Aqcd- Indeed, B.„ ^ 4GeV at our values of /3, 
so that aB^/2 ^ 1 at /5 = 6! Thus, although the r.h.s. of Eq. (p^) is formally of O(a^), it 
can be comparable in magnitude to an 0{a) effect. Of course, the numerator also depends 
on ca, although in a way which cannot be estimated simply. Thus, it is possible that the 
enhanced ca dependence cancels in Zy, and our results indicate that this is what happens: 
2- and 3-point discretizations lead to consistent results. This is an example of our general 
observation (see Sec. |V|) that the value of ca determined from Eq. ( jTSD improves the axial 
current in other correlation functions. 

In contrast, the improvement of the axial current does not guarantee that there are no 
enhanced 0(a) errors in the slope, by — bA [§]• In particular, using a mass dependent ca in 
{Aj)\ produces an enhanced higher order effect proportional to B.^. We see this clearly 
in our results. For example, for the 62NP data set, by — bA is —0.11(3) [—0.07(3)] for the 
chirally extrapolated ca and 2-point [3-point] discretization, while using the mass-dependent 
Ca these results change to —0.30(4) [+0.34(5)]. It is reassuring that the discretization de- 
pendence is much weaker for chirally extrapolated ca, since this is the choice we have made 
at the order of improvement that we are working (see Section ^. These are the results we 
quote. 
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FIG. 12: A fit of the form cy 



-V 



+ c^y /{mi - ma) to the 62NP data. 



VIII. Z\ 



The AWI which yields the best signal for Z\ is 



:(i2) 



^23) 



(31), 



Y.,{5sr {A,)r\y,y,)vrm 



\(13)/- 



T.A{Vir\y.y^)vrm 



(31), 



Z° (1 + 6vam3/2) 
Zl-Zl{l + hAam^/2) 



(25) 



which holds after extrapolation to rhi = m2 = 0. The intermediate state in these correlators 
is the vector meson. The quality of the signal for the ratio on the l.h.s. is illustrated in 
Fig. |T^. An example of the fit versus ^3/2 is shown in Fig. |TB|. The resulting values for 
Zy/lZ^)"^ and 6^ — by are given in Tabs. p| - |V^ . The latter have much larger errors than 



those in the determinations described in the previous section. 

Rather than obtaining Z^ by combining the results for Zy/lZj^)"^ with those obtained 
previously for Zy, it turns out to be better to use the product of the left hand sides of 
Eqs. (^Of ) and (p5|), which yields 1/(Z°)^ directly. Note that the hnear 0{a) fh^ dependence 
cancels in this product. The data, illustrated in Fig. |1^, show a dependence on m^ at the 2a 
level. Our final results are obtained by a fitting a constant to the data. A linear fit reduces 
the value by 1 — 2cr. 

As shown in Tabs. pI| - [VIIl| , the statistical errors in Z^ are roughly an order of magnitude 



larger than those in Zy, and are comparable to the size of the expected O(a^) terms. Thus, 
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FIG. 13: A constant fit as a function of m^/2 to extract cy from the 62NP data. Points included 
in the fit are superimposed with a fancy cross. 

either the statistical or the 0{a?) corrections can explain the difference between our results 
and those of the ALPHA collaboration, which are at the 1 — 2cr level. The deviations from 
1-loop perturbation theory are of the size expected if the 2- loop terms are ~ a^. 



IX. Z%/Z% AND bs - bp 

Our best estimates of Zp/{ZgZ^) and bs — bp are obtained from 






Z° (1 + bpam^/2) 



(26) 



with J = P 01 A^. Both numerator and denominator have pions as intermediate states, 



and have very good signals. Examples of their ratio are shown in Fig. |1^. As discussed 
in Sec. 0, this ratio should be independent of y^ up to higher order discretization errors. 
These errors are expected to be larger for the 60TI data set than for those with the non- 
perturbatively improved action, since the former are of Oia)^ and the latter of 0(a^\ Our 
results are qualitatively consistent with these expectations, as illustrated in Fig. |18|. Note 
that the scale is much finer for the lower graphs. A linear fit to Eq. (^61) gives our estimates 
for Zp/(Z^Z^) and bp — bs quoted in the tables. 
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FIG. 14: Linear fits to the l.h.s. of Eq. (20) versus (i) the AWI quark mass m (crosses) and (ii) 
the VWI quark mass m (diamonds). The data set is 62NP. The fit to crosses gives 6^ — by, while 
that to diamonds gives bA — by- 



Another way to extract Zp/(Z^Z^) is to use the relation between the two definitions of 



quark mass |]ll 



m 



m 



70 70 



{hA - bp)amav + b^a 



,m 



m„ 



(27) 



where Xav = (-^i+Ar2)/2. This relation is useful because Z^ = l/Z'^ and bs = —2bm |^, |24 |. 
From the non-leading terms one can, using non-degenerate quarks, separately determine 
bA — bp and bm- In this section we discuss and use only degenerate quarks, from which one 
can determine bA — bp — bm- The use of non-degenerate quarks, which allows a separate 
determination of 6^ ~ bp and bm, is discussed in Sec. PO. 



We have analyzed Eq. ( P?] ) by extracting fh from Eq. (^) using both the mass dependent 
and chirally extrapolated values of ca- An example of the data and linear fits is shown in 
Fig. |TP|. The intercepts are consistent, and we quote, in Tabs. |TT1|-|VI|, the results using the 
mass-dependent ca- We also show in the same figure the fit to Eq. (PB|), which should have 
the same intercept up to 0{a^) terms. While the data show no significant discrepancy, 
the results from Eq. (^) can have enhanced discretization errors. Indeed, it follows from 
Eq. (|15D that a change Aca results in 



^^ = AcAa-—:r = AcAaB^, 
m 2m 



(28) 
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FIG. 15: Illustration of the signal for the ratio defined in Eq. (|2 
propagators in all cases. 



for the four data sets, using ^3 



This is also the fractional change in the result for Zp/{ZgZ^) obtained from Eq. (p7|). Since, 
as noted above, 5^ ~ 4 GeV, this nominally 0{a^) uncertainty can be enhanced. As a result, 
we consider the evaluation using Eq. (pTl) less reliable than that based on Eq. (|^), and we 
use the latter as our best estimate. _ _ _ 

The slope of the linear fits to Eq. (|27|) for degenerate quarks gives Ba — bp + bs/2. The 
statistical errors on the results are small, but there is a systematic dependence on whether 
we use the mass- dependent or chirally extrapolated ca — an 0{a) effect enhanced by B.„. 
This problem is clear from Fig. |1^, and to highlight the magnitude we quote both values in 
Tabs. p^[VI| : the first corresponds to the mass- dependent ca and the second to the chirally 
extrapolated ca- Unlike the case of Ba — by, here the mass-independent ca, which is our 
choice, leads to results which depend very strongly on the choice of discretization. Because 
of these very large O(a^) effects, we do not use these estimates any further. 



Our derived results for Zp/Zg, presented in Tab. Vlll , are significantly smaller than the 



predictions of 1-loop perturbation theory. As noted in Sec. |IV|, the difference can only be 
explained by an unlikely 2-loop contribution ~ 4a^. 
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FIG. 16: Linear fit to the ratio defined in Eq. ( pq ) after extrapolation to mi = m,2 
62NP data set. 



for the 



X. ct 



To determine ct we consider the AWI for the bilinear Tij, i.e. 



(29) 



As was the case for cy, tuning Ct in order to make the ratio independent of y^ does not 
work. Instead we rewrite the identity in the following form, 



(12) ^(23), 



.(31), 



1 + acT 



EyA[-dM''Hy,yMrm .Eyi^sr W^^y^ ^^(o)) 



.(13) 



(31)/ 



Y.y{ni^\y.y^nTm 



^^ j:y{Tif\y,y,)Tif\o)) 



(30) 



where we have moved the ct dependence in (Tj)m onto the l.h.s., and used the fact that 
{Ti)ij has no contribution from the Ct term at p = 0. Given Z^, Eq. (^) determines Ct 
after the mi -^ extrapolation. The data for the ratios on the left and right hand sides 
of Eq. ( ^OD are illustrated in Figs. ^ and ^ respectively and expose the reason for the 
failure to extract ct by tuning with respect to 7/4: the two ratios are essentially flat within 
the domain of the chiral rotation (which roughly corresponds to the region of the fits in the 
Figure) . 
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FIG. 17: Z^, obtained from the product of ratios of correlators defined in the l.h.s. of Eq. (|20| ) and 

K5 points, as indicated by 



Eq. ([25|), shown as a function of m^/2. The constant fit is to the K2 
the fancy crosses. The data set is 62NP. 



ct should be independent of ms, up to corrections of 0{a^). Our results are consistent 
with this expectation at the 1 — 2cr level, as illustrated in Fig. ^ for the 62NP data set. 
Our quoted results are the weighted average over the ^2 — ^^5 points. 

To extract hx using the method proposed in [Q requires studying this AWI with all three 
quarks in Eq. (|30|) having different masses. We have not done this extended calculation, and 
consequently have no results for bx- 



XI. ADDITIONAL RELATIONS 



There are two additional relations that can be used to obtain information on improvement 

The first is 



constants. These were derived in Ref. 11, and discussed further in Ref. 



hp -hA 



4mi2 - 2 [mil +"^22] 
a[mii - m22]^ 



(31) 



An illustration of our results for the r.h.s. is shown in Fig. ^ and the results from fits to a 
constant are collected in Tabs. |IT^|VI|, and used to obtain the final results given in Tab. |V111 . 
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FIG. 18: Comparison of the signal in the ratio of correlators on the l.h.s. of Eq. 
Zp/Z^. The data are for k^ propagators in all cases. 



used to extract 



The second relation is 



bs-bv 



+ {bp - b, 



Ai2 - Rz[mii -m22] 



ni 



22J 



Aio = 






R, 






yO 



yO yO 

Zjp Zjy 



(32) 
(33) 
(34) 



As discussed in Ref. |], of the two kinds of sources: J^^^^ = ^^P(^^)(2', Z4)P(^^)(0), and 
j(2i) _ ^(21) ^[f^i^ < t < Zi that one can use in Eq. (^5]), the first has a better signal and 
smaller discretization errors. Unfortunately, the final results, quoted in the Tables, have 
very large errors due to large cancellations between the terms in the numerator on the r.h.s.. 
We, therefore, do not use this second combination in our final extraction of the individual 
6's given in Tab. |Vlli. 
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FIG. 19: Comparison of the quality of the Hnear fits used to extract Z^Z^/Zp. The three fits 
correspond to (i) Eq. (p6|), (ii) Eq. (^) with fh defined using the mass dependent ca, and (iii) 
Eq. (27) with fh defined using the chirahy extrapolated ca_- The data are from the 62NP set. 
Note that the intercepts from all three fits should agree up to errors of 0{a'^), but the slope of (i) 
is bp — bs whereas those of (ii) and (iii) are bp — bA — bs/2. 



XII. EQUATION-OF-MOTION OPERATORS 

The method for calculating the combination c'p + c'q of coefficients of equation-of-motion 
operators has been described in Sec. 0. The calculation, using Eq. (^^, involves three pieces. 
The slopes so are obtained from a linear fit to the l.h.s. of Eq. (|1^) versus fhi at fixed fh^. 
Examples of these fits are shown in Fig. ^, for the 62NP data set. The on-shell quantities 
Xo{bso — bo) and XobA can be obtained by combining results discussed in previous sections. 
The results for these three contributions, for the 62NP data set, are collected in Tab. ^. 
We find that XoBa gives almost the entire contribution. The final estimates for individual 
equation-of-motion constants are given in Tab. [TX . 

We briefiy discuss some details of the calculation, and the quality of the signal, in each 
of the five cases. 

• c'p + c'y\ We choose J = P and O = V4 {60 = A4), in which case the intermediate 



state is a pseudoscalar. 
c'p + c'j^: We choose J = 



Vi and O = Ai {60 = Vi) whereby the intermediate state is 
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FIG. 20: The signal in the ratio of correlators defined on the left hand side of Eq. (^) using K3 in 
all quark propagators. 

a vector meson. Unfortunately, the uncertainty in cy feeds in through 60 = Vi and 
affects the extraction of sa- Thus, even though sa contributes little to the central 
value of dp + c^, as illustrated in Tab. ^, it dominates the error. 

• 2cp: We choose J = S and O = P {50 = S). In this case, the intermediate state is a 
scalar and the signal is poor.^ The largest part of the error in dp comes from sp. The 
resulting uncertainty in dp dominates the error in the final estimate of dy, c^, and c^. 

• dp + dg-. The choice J = P and O = S {60 = P) gives a good signal in the correlation 
functions as the intermediate state is pseudoscalar. 

• dp + drp-. We choose J = T/.^ and O = Tij {60 = T^^). All correlation functions have 
a good signal as the intermediate state is a vector meson. 

The signal for sy, s^, and st is good for all ms, and leads to a reliable estimates with 
comparable errors for dp + Cy, dp + c'^, and dp + c^. In all cases we find that sq are 



^ A better choice might be to use J — 'Yl,sP{'^)P{^) with Z4 ^> 2/4 3> 0, in which case the intermediate state 
is pseudoscalar, but this requires an extra inversion. 
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FIG. 21: The signal in the ratio of correlators defined on the right hand side of Eq. 
in all quark propagators. 



using ^3 



c'a + C'p 


so 


-^^0(^50 - ^o)/2 


XobA 


Cy + Cp 


-0.27(04) 


-0.07(2) 


1.52(4) 


C'^ + C'p 


-0.13(06) 


0.07(2) 


1.41(3) 


Cp -\- Cp 


-0.73(16) 


0.06(2) 


1.70(7) 


c's + c'p 


-0.14(03) 


-0.05(1) 


1.29(3) 


c^ + c'p 


-0.23(05) 


0.02(3) 


1.46(4) 



TABLE X: The three contributions to the coefficient of the equation of motion operators c'q + c'p 
for the 62NP data set. 

independent of m^ within statistical errors. Our final results are given by the weighted 
mean over 777,3 corresponding to K2 — K5. 

To compare to the predictions of perturbation theory, it is best to use the results for 
c'x + c'p, X = V,A, S, T, in the upper part of Table |^ since these have the smallest 
statistical errors. These four quantities are indeed consistent with the expected result 2[1 + 
0{as) + 0{a)]. The fifth quantity, 2cp, is only determined reliably at /? = 6.2, and also 
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FIG. 22: Estimates of ct as a function of raj^jl for the 62NP data set. The constant fit is to the 
1^2 — n^ points. 

agrees with this expectation. These agreements are a consistency check on the extension of 
the improvement program to off-shell quantities. 

XIII. CONCLUSION 

We have demonstrated the feasibility of the WI method, with non-degenerate quark 
masses, for determining the improvement and scheme-independent normalization constants 
of the quark bilinear operators. The main advantage of using non-degenerate quarks is that 
one can extract all the hx- These quantities effect the overall normalization of operators away 
from the chiral limit, and their determination is relevant to phenomenological applications 
involving heavy mesons. 

Our implementation of the Ward identities differs substantially from that used by the 
ALPHA collaboration, so that the results from the two methods can differ. These differences, 
should, however, be of size 0{a) and 0{a?), respectively, for improvement and normalization 
constants. The differences between the two sets of results are, in fact, consistent with these 
expectations. We stress, however, that for the small quantities, ca and cy, this "consistency" 
allows a substantial uncertainty at /? = 6. For example, Ac^ = 0.05 would lead to an k. 10% 
uncertainty in /^r and 3% in /^i. At /3 = 6.2, on the other hand, there is a much smaller 
variability. 
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FIG. 23: A constant fit to the 62NP data for bp - Ia obtained using Eq. (|3]j 



Both cy and ct are obtained as a small difference between two large terms. We are, 
nevertheless, able to extract these quantities with reasonable precision. In particular, in the 
case of cy, we find that our best results come from enforcing a different Ward identity than 
considered previously, with a consequent reduction in errors. This improvement is important 
for phenomenological applications (see, e.g., Ref. ^), and also leads to smaller errors in our 
results for Z% Zp/Zg, ct and c^. 

On the whole, tadpole-improved 1-loop perturbation theory underestimates the deviations 
of renormalization and improvement constants from their tree level values. In all but one 
case, however, these discrepancies can be understood as a combination of a 2-loop correction 
of size (1— 2)xa^ [for Zy, Z% andcAJ, higher order discretization errors of size (1— 2)xaAQCD 
[for cy, Ct and by], and statistical errors [for 6^, bp, and 65]. The only exception is Zp/Zg, 
for which a very large higher order perturbative contribution of size 4 x a^ is needed to 
reconcile our non-perturbative results with 1-loop perturbation theory. 

We have, for the first time, presented results for the coefficients of equation of motion 
operators that are needed to improve the theory off-shell. The most striking feature of their 
calculation is the improvement in the reliability of the calculation between f3 = 6.0 and 6.2. 

An important issue is at what quark mass 0{a) improvement breaks down, due to our 
neglect of higher order terms. To address this issue we examine the case of the charm quark 
at /9 = 6.2 for which ma ^ 0.5 and rh ^ 0.4. Since bx ~ 1-1, the 0{a) corrections to Zx are 
approximately 45%. Assuming geometric growth, this would imply ^ 20% correction from 
the neglected Ola"^) terms. This is indeed what we find for Zy, for which non-perturbative 
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FIG. 24: Linear fits to the l.h.s. of Eq. (10), the slopes of which, sq, determine the coefficients of 
the equation of motion operators. The data set is 62NP and fh^ corresponds to K3. 



results for charm quarks are available, and the data are good enough to allow the quadratic 
fit given in Eq. (0). On the other hand, we find that if we use the alternative 0{a) improved 
expression Zy = Zy{l + byma), it works to within 1% at the charm quark mass. 

Finally, we stress that the use of non-degenerate quarks to determine the bx and ct 
could be applied equally well in the context of the Schrodinger functional. It would be very 
interesting to compare results so obtained to those we have found here. 
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APPENDIX A 

In this appendix we review the relation between continuum and lattice fields and the 
1-loop perturbative results. Throughout this paper we use 

(^/?) continuum = v4«^(Ci?)lIttice • (Al) 

This normalization makes comparison between tadpole-improved 1-loop and non- 
perturbative results, quoted in Tab. |V111| , straightforward. In the tadpole improved the- 



ory [^, the normalization commonly used is -^AkiKjU^. To maintain the field normalization 
as ^JAKiKj we have absorbed uq into Zq^^^.^. Consequently, the TI perturbative result we 
use is Zqp^j.^ = Mo(l + to(y^^), where to is the TI 1-loop coefficient. 
A second way in which tadpole improvement is defined is 



^fo) _ Q^ . /i ^ ^. /i ^ ^^0 rrn\iii) 



(C^/j) continuum " S^cyl 'aZ~ V ^ 'T~ ^ O ,pert{^) lattice 



where we have again absorbed a factor of Uq in Zq ^ to maintain the same definition as 
above. Eq. ( |A^ ) shows that using tadpole-improved field renormalization is equivalent, at 
0{a), to using ho = Skc in Eq. (||). In tree-level TI perturbation theory 8kc = ^/uq, and is 
the appropriate value for bo as shown in Eq. (|A6|) . 

The 1-loop perturbative calculations have been done by the ALPHA and JLQCD col- 
laborations [^^ ^, ^. Here we express the results for the tadpole improvement scheme 
stated above. Tadpole improvement requires choosing a quantity, Mq, which is unity at tree- 
level, whose perturbative series is dominated by a tadpole contribution, and which can be 
evaluated non-perturbatively. Any other quantity X, whose perturbative expansion is 

X = X(°)+X«a,, (A3) 

and which is dominated by n contributions of the tadpole diagram, can then be re-written 



as 



Xti = <(X(o)+x(!;a,,T/), 
where 

X^'} = X(i) - nX^^^u'i^ . (A4) 

Here Mq is the coefficient of a^ in the perturbative expansion of Mq, and ag^ri is an improved 
coupling that we choose to be g'^/iiruQ, where j3 = G/g"^. Since all results we quote are 
tadpole- improved, we henceforth omit the subscript TI for brevity. 

In this paper, we choose, for uq, the fourth root of the expectation value of the plaquette 
for which Mq = — vr/S. Our Monte Carlo data yields uq = 0.8778 at /3 = 6.0 and 0.8851 at 
P = 6.2. Using this uq, we find that as = 0.1340 and 0.1255 at the two /3's. 
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7r 




Cp 


Op 


I? 


s 


1 


-1.002 




1.3722 


1.2818 


p 


1 


-1.328 




0.8763 


0.7859 


V 





-0.579 


-0.2054 


0.8796 


0.7892 


A 





-0.416 


-0.0952 


0.8646 


0.7742 


T 


-4/3 


-0.134 


-0.1505 


0.7020 


0.6116 



TABLE XI: The tadpole-improved one-loop coefficients in Eq. (A6). The tadpole- improvement 
factor, uq, has been chosen to be the fourth root of the plaquette expectation value. 



At one- loop, the coefficient of the clover term is 

CSW = M(7^(l + c'swOis) 



(A5) 



where Csl/ = 0.214 is obtained by converting the results by Wohlert |3^ and the ALPHA 



-sw 



collaboration |^] to tadpole-improved form. Then csw = 1.521 and 1.481 at /5 = 6.0 and 
6.2 respectively. 

The tadpole-improved renormalization constants at one loop are given by the formulae: 



,7r 



.(!)> 



Z° = uo[l + as{j^HH' + 4")] 
47r 



cr 



a.oC 



(1) 



s'^r 



6r = [1 + a^?] 



(A6) 



where /i is the scale at which the continuum MS theory is defined. The final results for all 
these tadpole- improved coefficients are given in Tab. ^. There are two points worth noting: 
(i) the tadpole factors cancel in the product b-pm, whereas neither br nor rh has any; (ii) 
the one- loop correction, Cg^y, does not contribute to the renormalization or improvement 
constants at 0{as). 

APPENDIX B 



In this appendix we review tree-level improvement of Wilson fermions and define our 
conventions for improvement coefficients. The 0{a) improvement of Wilson fermions can be 



obtained by the transformation [^ 






(IT ^ 



^ 






(Bl) 
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where the continuum equation of motion is given by {p + m)ip = 0. Using the fact that the 
Wilson-clover operator aW is related to p by 

aWtp = a(j/) + m)ij + 0{a^) 

V^aVV = HP - m)a + 0{a^) , (B2) 

we can rewrite the improved fermion fields ipi and ipi as 

?/ = V^ {l + ^[cs^rP + (2 - c,.^r)m] + «^(1-^^^^) ^| + o(a2) , (B3) 

where Cswr represents an arbitrary 'rotation' parameter. Operators composed of these im- 
proved fermion fields are automatically 0{a) improved at tree level. 

In particular, we can construct the tree-level improved fermion bilinears S, P, V, A and 

T, as 

Si = {l+armbs)SL-arcsSi.nnk J^^s 



arc' 



Pi = (1 + armhp)PL + arcpdi^Ai^t, 'I~^p 



arc' 



Vi,^ = (1 + armbv)VL,^, + arcydyTL^^^y - arcy \/i_iink,/. -r—E- 



V,ti 



Ai,f, = (1 + armbA)A^^ + arcAd^^PL - arcAAi.n^],^^, J^^^-'^ 



arc' 



Ti,f,u = (1 + armbT)TL,f,u + arcTid^Vi^y - d^^V^f,) - arcrTi_Hnk,^^, ^^t,mi^ , (B4) 

where, we have dropped all 0{a'^) terms, and for all O, bo = {2 — Cs^r)/2, Cq = Cgwrl'^ 
(except ct = —Cswr/'^)^ c'^ = 1 — Cgwr and co = Cgwr/'^- The local operators, Ol, are defined 
as ^To'ip with To being 1, 75 = 71727374, 7^' 7^75 and ia^,^ = -[7^, 7^]/2 for O = S, P, V^, 

A^ and T^,^ respectively^; the equation of motion operators, Eo, as ipiXoW ~ W^o)"^, and 
the 1-link operators Oi-imk as 

5'l-link = i'pip 
^l-link,/. = V'-D/i^ 
Al.\ink,f, = -i'ipDuC^ufj.lb'^ 

Ti.imk,f,u = e^yX5'4'Dx'j575'^, (B5) 

where D = D — D- It is easy to see that the operators O/,, Oi.Hnk, and Oem form an over- 
complete basis for all dimension-4 fermion bilinear operators, and therefore no new operators 



In Ref. H, a factor of i was inadvertantly missed in the definition of (7^,y. The correct definition is 
CT^i. = i[7^,7,.]/2. 

42 



are needed for non-perturbative improvement of the quenched theory. In this paper, we have 
chosen to ehminate the 1-hnk operators (and the <9^A^ term in Pj) non-perturbatively by 
an appropriate choice of Cswr- At tree-level, this implies Cswr = 0, whereby 



bo = 


= 1, 


Co = 


= 0, 


c'o- 


= 1. 



(B6) 

It is important to note that beyond tree-level, the matrix elements of the 1-link operators 
have divergences proportional to a~^, and hence contribute to the renormalization constants 
at 0{a'^). As a result, not only do the 0{a) correction terms bo, co and c'q depend on the 
choice of Co, but so do Z^, except for O = P. 



APPENDIX C 

In this appendix we give a brief description of the two exceptional configurations we 
found in the 60NP data set. In both of these we find that the zero mode is localized over 
5 — 10 timeslices. If the Wuppertal source overlaps with the zero mode then the norm of 
the pion propagator for quark mass kj can be up to a factor of a hundred larger than the 
average over the remaining configurations. If, on the other hand, the source time slice does 
not overlap with the zero mode, then we observe a "normal" temporal fall-off in the pion 
correlator until it hits the zero mode, when it shows a large bump. These two anomalous 
behaviors are illustrated in Fig. |2^. 
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and the large deviation from exponential fall-off if it does not. In each case the time coordinates 
are translated so that the Wuppertal source is at t = 1. 
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